# This file and accompanying functions are renamed with domestic RC added to 
# a CF on BLP data
#
rm(list=ls())
#library(data.table)
library(HeadR)
#devtools::install_github("ckhead/HeadR")
#remotes::install_github("ckhead/HeadR")
#library(doParallel)
#library(doRNG)
library(SQUAREM)
library(fixest)
library(tictoc)
#library(xtable)
#library(latex2exp)
setwd("~/Dropbox/BLP_new")
source("Rfiles/CF_BLPdata_domestic_functions.R")
source("Rfiles/CF_MMX_functions.R")
source("Rfiles/CF_MMX_MLOG_functions.R")
library(Rcpp)
#check that Rcpp works on this computer
evalCpp("2 + 2")
# create the Rcpp function we need for computing equilibrium
Rcpp::sourceCpp("Cppfiles/crossDelta.cpp")
#parallel processing set up for doParallel, doRNG
#cl <- detectCores()-1
#registerDoParallel(cores=cl)
#getDoParWorkers()
MAXIT <- 500
seedn <- 777 # 666
sequences <- "normal" # alternatives: normal, halton or "sobol"
#
n.hh <- 100000 # careful, this is slow, but works
#n.hh <- 1000 # fast & dirty
U0 <- 1  # needed in prob.mat(), but not prob.mat.namedxvars()
n.reps <- 1 # for richsubs (no need to do a lot if a lot of consumers.)
sigmadomesticpct.list <- c(0,3,25,34,50,72,100,125,150,172,200,225,250) # either  0, 72 or 172
num.loops <- length(sigmadomesticpct.list)
#
# This is intended to choose which country is targeted
source(file = "Rfiles/CF_BLPdata_brandiso.R")
target.list <- brands[iso!="USA",unique(iso)] # all non US brands
#target.list <- "JPN" # the countries to which the tariff is applied.
rm(brands)
#
#
# Here choose which set of parameters to use. 
param.all <- fread("Data/BLP99/pyblp_domestic_param.csv") # estimated using C&G package
param.list <- names(param.all[,par_abbrev:=NULL])
param.list 
#add pub_param (the BLP 99 parameters)
#param.list <- c("pub_param","Domestic_pub_param",param.list)
param.chosen <- c("Domestic_pub_param")
#
#
tic()
#Loop over sigmadomesticpct list, conducting counterfactuals ====
for(sigmadomesticpct in sigmadomesticpct.list[4:num.loops] ) {
cat("Now running for sigma domestic =", sigmadomesticpct ,"\n")
#
# This is used in the counterfactuals : number of characteristics (outside the constant)
if(param.chosen %like% "^Domestic"){
n.x <- 5 #THIS NEEDS TO BE 5 WHEN ADDING DOMESTIC, 4 otherwise (we could automate)
} else n.x <- 4
#
###
# What is the average own elasticity? ====
###
setwd("~/Dropbox/BLP_new")
source("Rfiles/CF_BLPdata_domestic_getdata.R") # get original data: original published_param.csv is 1995 and should not be changed (used in other files)
#
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#initialization phase
set.seed(seedn)
data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)] # super assignment needed here because of scope for DGP.BLP
BLP.meanownelas
#
###
# Implement CF on BLP ====
###
#
source("Rfiles/CF_BLPdata_domestic_getdata.R") # get original data
year.gph <- 1990
tar <- 0 # original tariff (to be raised)
tar.increase <- 0.1 # the actual increase in tariff rate
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix 
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#do a settings file only for the speed 
settings.var <- data.frame(speed =c(0.3,0.3,0.3))
#
# Run the calibration to get xi and mc
year.monte <- year.gph
objective <- "CF"
market <- 1
#
#Run sims for the 3 versions
#EHA:
set.seed(seedn)
CF_BLPdata.simulv2 <- CF_BLPdata.simul(vsim=2,cf.meth="EHA")
set.seed(seedn)
CF_BLPdata.simulv3 <- CF_BLPdata.simul(vsim=3,cf.meth="EHA")
#
CF_BLPdata.simulallv <- rbind(CF_BLPdata.simulv2,CF_BLPdata.simulv3) 
CF_BLPdata.table <- CF_BLPdata.simulallv[,list(sigmadomesticpct,v, chg_true_qty,chg_ces_qty)]
#
saveRDS(CF_BLPdata.table,paste0("Tables/CF_BLPdata/SD/CF_BLPdata_simul_",param.chosen,sigmadomesticpct,".rds"))
system('echo "Done with one loop" | mail -s "updatefrom Gonzalez" ckhead@gmail.com')
}
toc()

#


